
## make sure to run the codes in pica2_pure.Rmd before running these codes 

rm(list = ls())
library('stringr')
library(ggplot2)
library(tidyverse)

load("pica2_political_final.rdata")
load("pica2_econ_final.rdata")

## #################
## POLITICAL Models 
## #################
pol_list <- political_list
PA_ci_l_us <- pol_list$Effect_US[c(2,3), 5] - pol_list$Effect_US[c(2,3), 4]
Bound_ci_l_us <- pol_list$Effect_US_bound[c(1,2), 4] - pol_list$Effect_US_bound[c(1,2), 3]
PA_ci_l_china <- pol_list$Effect_China[c(2,3), 5] - pol_list$Effect_China[c(2,3), 4]
Bound_ci_l_china <- pol_list$Effect_China_bound[c(1,2), 4] - pol_list$Effect_China_bound[c(1,2), 3]

mean(c(PA_ci_l_us, PA_ci_l_china)/c(Bound_ci_l_us, Bound_ci_l_china))

Effect_US_pol <- rbind(pol_list$Effect_US[c(1,2,3), c(1, NA, NA, 4, 5)], 
                       setNames(cbind(NA, pol_list$Effect_US_bound[1:2, c(1:4)]), 
                                colnames(pol_list$Effect_US[c(1,2,3), c(1, NA, NA, 4, 5)])))
Effect_China_pol <- rbind(pol_list$Effect_China[c(1,2,3), c(1, NA, NA, 4, 5)], 
                          setNames(cbind(NA, pol_list$Effect_China_bound[1:2, c(1:4)]), 
                                   colnames(pol_list$Effect_China[c(1,2,3), c(1, NA, NA, 4, 5)])))

pol <- rbind(Effect_US_pol, Effect_China_pol)
pol <- as.data.frame(pol)
rownames(pol) <- NULL
colnames(pol) <- c("Estimate", "l_bound", "u_bound", "cilower", "ciupper")
rownames(pol) <- c("ATE_US", rownames(pol_list$Effect_US)[2:3], paste0("bound:", rownames(pol_list$Effect_US_bound)),
                   "ATE_China", rownames(pol_list$Effect_China)[2:3], paste0("bound:", rownames(pol_list$Effect_China_bound)))

pol$group <- c(rep("US", 5), rep("China", 5))

pol$row_names2 <- c("ATE", "Choosers", "non-choosers", "Choosers", "non-choosers",
                    "ATE", "Choosers", "non-choosers", "Choosers", "non-choosers")

pol$row_names2 <- factor(pol$row_names2, level = c("ATE", "Choosers", "non-choosers"), ordered=TRUE)
pol$row_names2 <- str_wrap(pol$row_names2, width = 11)
pol$group <- str_wrap(pol$group, width = 10) # adjust width as per your requirement
pol$type <- c('ATE', rep('ACTE', 4), 'ATE', rep('ACTE', 4))
pol$pica_style <- c('none', rep("pica_2", 2), rep("pica_1", 2),
                    'none', rep("pica_2", 2), rep("pica_1", 2))

dodge_width <- 0.5
pol$row_num <- case_when(pol$row_names2 == "ATE" ~ 1, 
                         pol$row_names2 == "Choosers" ~ 2,
                         pol$row_names2 == "non-choosers" ~ 3)

pol_US    <- pol[pol$group == "US", ]
pol_China <- pol[pol$group == "China", ]

f_pol_China <- ggplot(data = pol_China, aes(x = row_num, y = Estimate, shape = pica_style, colour = pica_style)) +
  geom_point(data = pol_China[1,], aes(x = row_num, y = Estimate, shape = pica_style), size = 4, # shape = c(16), 
             position = "identity") + 
  geom_errorbar(data = pol_China[1, ], aes(x = row_num, ymin = cilower, ymax = ciupper), width=.1, linewidth = 0.5,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  geom_point(data = pol_China[2:3, ], aes(x = c(row_num - 0.1), y = Estimate, shape = pica_style), size = 4, # shape = c(17), 
             position = "identity") + 
  geom_errorbar(data = pol_China[2:3, ], aes(x = row_num - 0.1, ymin = cilower, ymax = ciupper), 
                width=.1, linewidth = 0.5,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  geom_errorbar(data = pol_China[4:5, ], aes(x = row_num + 0.1, ymin = l_bound, ymax = u_bound), width=.13, linewidth = 2,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  geom_errorbar(data = pol_China[4:5, ], aes(x = row_num + 0.1, ymin = cilower, ymax = ciupper), width=.1, linewidth = 0.5,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  scale_x_discrete(limits = c("ATE", "Choosers", "non-choosers"),
                   labels = c('ATE', 'ACTE:\n Choosers', 'ACTE:\n Non-Choosers')) + 
  geom_hline(yintercept = 0, linetype = 3) +
  scale_shape_manual(labels = c("ATE", "ACTE (PICA-2)", "ACTE (PICA)"), values=c(16, 17, 15), 
                     breaks = c("none", "pica_2", "pica_1"),
                     guide = guide_legend(override.aes = list(linetype = c(1, 1, 1),
                                                              shape = c(16, 17, NA),
                                                              color = c("black", "red", "blue"),
                                                              linewidth = c(1, 1, 2)))) + 
  scale_colour_manual(labels = c("ATE", "ACTE (PICA-2)", "ACTE (PICA)"), 
                      values=c("black", "red", "blue"), breaks = c("none", "pica_2", "pica_1")) + 
  facet_wrap(. ~ group, 
             # scales = "free", space = "free", 
             ncol = 1,
             labeller = as_labeller(c("China" = "Effects of Chinese Propaganda", "US" = "Effects of US Propaganda"))) +
  theme_bw() + 
  theme(# legend.title = element_blank(),
    # legend.text=element_text(size = 16),
    # legend.position = c(0.25, 0.83), 
    legend.position = "none",
    legend.title = element_blank(),
    legend.text=element_text(size = 14),
    axis.title.y = element_text(color = 'black', size = 14),
    axis.text.y = element_text(color = 'black', size = 12),
    axis.text.x = element_text(angle = 0, hjust = 0.5, vjust= 1, color= 'black', size =14, face ="bold" ), 
    panel.border = element_rect(colour = "black", fill= 'NA', size=1),
    strip.background = element_rect(colour = "black", fill = 'grey89', size = 1),
    strip.text = element_text(face = "bold", size = 13),
    plot.title = element_text(face = "bold", size = 14))  +  # adjust text properties here
  labs(title = "(A) Outcome: Prefer China's Political Model to US's",
       x = "",
       y = "Estimate"
  ) + ylim(c(-1.5, 3.5)) 

f_pol_US <- ggplot(data = pol_US, aes(x = row_num, y = Estimate, shape = pica_style, colour = pica_style)) +
  geom_point(data = pol_US[1,], aes(x = row_num, y = Estimate, shape = pica_style), size = 4, # shape = c(16), 
             position = "identity") + 
  geom_errorbar(data = pol_US[1, ], aes(x = row_num, ymin = cilower, ymax = ciupper), width=.1, linewidth = 0.5,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  geom_point(data = pol_US[2:3, ], aes(x = c(row_num - 0.1), y = Estimate, shape = pica_style), size = 4, # shape = c(17), 
             position = "identity") + 
  geom_errorbar(data = pol_US[2:3, ], aes(x = row_num - 0.1, ymin = cilower, ymax = ciupper), 
                width=.1, linewidth = 0.5,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  geom_errorbar(data = pol_US[4:5, ], aes(x = row_num + 0.1, ymin = l_bound, ymax = u_bound), width=.13, linewidth = 2,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  geom_errorbar(data = pol_US[4:5, ], aes(x = row_num + 0.1, ymin = cilower, ymax = ciupper), width=.1, linewidth = 0.5,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  scale_x_discrete(limits = c("ATE", "Choosers", "non-choosers"),
                   labels = c('ATE', 'ACTE:\n Choosers', 'ACTE:\n Non-Choosers')) + 
  geom_hline(yintercept = 0, linetype = 3) +
  scale_shape_manual(labels = c("ATE", "ACTE (PICA-2 Estimator)", "ACTE (PICA Bound)"), values=c(16, 17, 15), 
                     breaks = c("none", "pica_2", "pica_1"),
                     guide = guide_legend(override.aes = list(linetype = c(0, 0, 1),
                                                              shape = c(16, 17, NA),
                                                              color = c("black", "red", "blue"),
                                                              linewidth = c(1, 1, 2)))) + 
  scale_colour_manual(labels = c("ATE", "ACTE (PICA-2 Estimator)", "ACTE (PICA Bound)"), 
                      values=c("black", "red", "blue"), breaks = c("none", "pica_2", "pica_1")) +
  facet_wrap(. ~ group, 
             # scales = "free", space = "free", 
             ncol = 1,
             labeller = as_labeller(c("China" = "Effects of Chinese Propaganda", "US" = "Effects of US Propaganda"))) +
  theme_bw() + 
  theme(legend.position = c(0.65, 0.83), 
        legend.title = element_blank(),
        legend.text=element_text(size = 14),
        axis.title.y = element_text(color = 'black', size = 14),
        axis.text.y = element_text(color = 'black', size = 12),
        axis.text.x = element_text(angle = 0, hjust = 0.5, vjust= 1, color= 'black', size =14, face ="bold" ), 
        panel.border = element_rect(colour = "black", fill= 'NA', size=1),
        strip.background = element_rect(colour = "black", fill = 'grey89', size = 1),
        strip.text = element_text(face = "bold", size = 13),
        plot.title = element_text(face = "bold", size = 14))  +  # adjust text properties here
  labs(title = "",
       x = "",
       y = "Estimate"
  ) + ylim(c(-1.5, 3.5)) 

# ###################
# Economic Models 
# ###################
Effect_US_econ <- rbind(econ_list$Effect_US[c(1,2,3), c(1, NA, NA, 4, 5)], 
                        setNames(cbind(NA, econ_list$Effect_US_bound[1:2, c(1:4)]), 
                                 colnames(econ_list$Effect_US[c(1,2,3), c(1, NA, NA, 4, 5)])))
Effect_China_econ <- rbind(econ_list$Effect_China[c(1,2,3), c(1, NA, NA, 4, 5)], 
                           setNames(cbind(NA, econ_list$Effect_China_bound[1:2, c(1:4)]), 
                                    colnames(econ_list$Effect_China[c(1,2,3), c(1, NA, NA, 4, 5)])))

econ <- rbind(Effect_US_econ, Effect_China_econ)
econ <- as.data.frame(econ)
rownames(econ) <- NULL
colnames(econ) <- c("Estimate", "l_bound", "u_bound", "cilower", "ciupper")
rownames(econ) <- c("ATE_US", rownames(econ_list$Effect_US)[2:3], paste0("bound:", rownames(econ_list$Effect_US_bound)),
                    "ATE_China", rownames(econ_list$Effect_China)[2:3], paste0("bound:", rownames(econ_list$Effect_China_bound)))
econ$group <- c(rep("US", 5), rep("China", 5))
econ$row_names2 <- c("ATE", "Choosers", "non-choosers", "Choosers", "non-choosers",
                     "ATE", "Choosers", "non-choosers", "Choosers", "non-choosers")
econ$row_names2 <- factor(econ$row_names2, level = c("ATE", "Choosers", "non-choosers"), ordered=TRUE)
econ$row_names2 <- str_wrap(econ$row_names2, width = 11)
econ$group <- str_wrap(econ$group, width = 10) # adjust width as per your requirement
econ$type <- c('ATE', rep('ACTE', 4), 'ATE', rep('ACTE', 4))
econ$pica_style <- c('none', rep("pica_2", 2), rep("pica_1", 2),
                     'none', rep("pica_2", 2), rep("pica_1", 2))

dodge_width <- 0.5
econ$row_num <- case_when(econ$row_names2 == "ATE" ~ 1, 
                          econ$row_names2 == "Choosers" ~ 2,
                          econ$row_names2 == "non-choosers" ~ 3)

econ_US    <- econ[econ$group == "US", ]
econ_China <- econ[econ$group == "China", ]

f_econ_China <- ggplot(data = econ_China, aes(x = row_num, y = Estimate, shape = pica_style, colour = pica_style)) +
  geom_point(data = econ_China[1,], aes(x = row_num, y = Estimate, shape = pica_style), size = 4, # shape = c(16), 
             position = "identity") + 
  geom_errorbar(data = econ_China[1, ], aes(x = row_num, ymin = cilower, ymax = ciupper), width=.1, linewidth = 0.5,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  geom_point(data = econ_China[2:3, ], aes(x = c(row_num - 0.1), y = Estimate, shape = pica_style), size = 4, # shape = c(17), 
             position = "identity") + 
  geom_errorbar(data = econ_China[2:3, ], aes(x = row_num - 0.1, ymin = cilower, ymax = ciupper), 
                width=.1, linewidth = 0.5,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  geom_errorbar(data = econ_China[4:5, ], aes(x = row_num + 0.1, ymin = l_bound, ymax = u_bound), width=.13, linewidth = 2,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  geom_errorbar(data = econ_China[4:5, ], aes(x = row_num + 0.1, ymin = cilower, ymax = ciupper), width=.1, linewidth = 0.5,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  scale_x_discrete(limits = c("ATE", "Choosers", "non-choosers"),
                   labels = c('ATE', 'ACTE:\n Choosers', 'ACTE:\n Non-Choosers')) + 
  geom_hline(yintercept = 0, linetype = 3) +
  scale_shape_manual(labels = c("ATE", "ACTE (PICA-2)", "ACTE (PICA)"), values=c(16, 17, 15), breaks = c("none", "pica_2", "pica_1")) + 
  scale_colour_manual(labels = c("ATE", "ACTE (PICA-2)", "ACTE (PICA)"), 
                      values=c("black", "red", "blue"), breaks = c("none", "pica_2", "pica_1")) + 
  facet_wrap(. ~ group, 
             # scales = "free", space = "free", 
             ncol = 1,
             labeller = as_labeller(c("China" = "Effects of Chinese Propaganda", "US" = "Effects of US Propaganda"))) +
  theme_bw() + 
  theme(# legend.position = c(0.25, 0.83), 
        legend.position = "none", 
        legend.title = element_blank(),
        legend.text=element_text(size = 14),
        axis.title.y = element_text(color = 'black', size = 14),
        axis.text.y = element_text(color = 'black', size = 12),
        axis.text.x = element_text(angle = 0, hjust = 0.5, vjust= 1, color= 'black', size =14, face ="bold" ), 
        panel.border = element_rect(colour = "black", fill= 'NA', size=1),
        strip.background = element_rect(colour = "black", fill = 'grey89', size = 1),
        strip.text = element_text(face = "bold", size = 13),
        plot.title = element_text(face = "bold", size = 14))  +  # adjust text properties here
  labs(title = "(B) Outcome: Prefer China's Economic Model to US's",
       x = "",
       y = "Estimate"
  ) + ylim(c(-3.5, 3.5)) 

f_econ_US <- ggplot(data = econ_US, aes(x = row_num, y = Estimate, shape = pica_style, colour = pica_style)) +
  geom_point(data = econ_US[1,], aes(x = row_num, y = Estimate, shape = pica_style), size = 4, # shape = c(16), 
             position = "identity") + 
  geom_errorbar(data = econ_US[1, ], aes(x = row_num, ymin = cilower, ymax = ciupper), width=.1, linewidth = 0.5,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  geom_point(data = econ_US[2:3, ], aes(x = c(row_num - 0.1), y = Estimate, shape = pica_style), size = 4, # shape = c(17), 
             position = "identity") + 
  geom_errorbar(data = econ_US[2:3, ], aes(x = row_num - 0.1, ymin = cilower, ymax = ciupper), 
                width=.1, linewidth = 0.5,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  geom_errorbar(data = econ_US[4:5, ], aes(x = row_num + 0.1, ymin = l_bound, ymax = u_bound), width=.13, linewidth = 2,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  geom_errorbar(data = econ_US[4:5, ], aes(x = row_num + 0.1, ymin = cilower, ymax = ciupper), width=.1, linewidth = 0.5,
                position = position_dodge2(width = dodge_width), size = 1, stat = "identity") + 
  scale_x_discrete(limits = c("ATE", "Choosers", "non-choosers"),
                   labels = c('ATE', 'ACTE:\n Choosers', 'ACTE:\n Non-Choosers')) + 
  geom_hline(yintercept = 0, linetype = 3) +
  scale_shape_manual(labels = c("ATE", "ACTE (PA Estimator)", "ACTE (Bound)"), values=c(16, 17, NA), 
                     breaks = c("none", "pica_2", "pica_1"),
                     guide = guide_legend(override.aes = list(linetype = c(0, 0, 1),
                                                              shape = c(16, 17, NA),
                                                              color = c("black", "red", "blue"),
                                                              linewidth = c(1, 1, 2)))) + 
  scale_colour_manual(labels = c("ATE", "ACTE (PA Estimator)", "ACTE (Bound)"),
                      values=c("black", "red", "blue"), breaks = c("none", "pica_2", "pica_1")) +
  facet_wrap(. ~ group, 
             # scales = "free", space = "free", 
             ncol = 1,
             labeller = as_labeller(c("China" = "Effects of Chinese Propaganda", "US" = "Effects of US Propaganda"))) +
  theme_bw() + 
  theme(legend.position = "none", 
        axis.title.y = element_text(color = 'black', size = 14),
        axis.text.y = element_text(color = 'black', size = 12),
        axis.text.x = element_text(angle = 0, hjust = 0.5, vjust= 1, color= 'black', size =14, face ="bold" ), 
        panel.border = element_rect(colour = "black", fill= 'NA', size=1),
        strip.background = element_rect(colour = "black", fill = 'grey89', size = 1),
        strip.text = element_text(face = "bold", size = 13),
        plot.title = element_text(face = "bold", size = 14))  +  # adjust text properties here
  labs(title = "",
       x = "",
       y = "Estimate"
  ) + ylim(c(-3.5, 3.5)) 


library(gridExtra)
library(grid)

pdf("figures/Figure_2.pdf", height = 9, width = 12)
grid.arrange(grobs = list(f_pol_China, f_pol_US, 
                          f_econ_China, f_econ_US), ncol = 2)
dev.off()